home *** CD-ROM | disk | FTP | other *** search
/ Meeting Pearls 1 / Meeting Pearls Vol 1 (1994).iso / installed_progs / gfx / lise2.1 / lise / src / fpeak.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-03-31  |  9.1 KB  |  374 lines
  1. /* 
  2.    Find peaks:
  3. */
  4.  
  5. #include <stdio.h>
  6. #include <spec.h>
  7. #define TMPFILE "global.def"
  8. #define MAXPEAKS 30
  9.  
  10. float *spc,*err,*tim;
  11. int   tri=0,
  12.       range = 1024;
  13.  
  14. help()
  15. {
  16.   printf("finding prompt peaks in file (first approximation)\n");
  17.   printf("fprompt file [options]\n");
  18.   printf("print out prompt peak positions\n");
  19.   printf("options:\n");
  20.   printf("  -tri [1/2] triangular spectrum beginning with minimum: 1\n");
  21.   printf("                                                peak:    2\n");
  22.   printf("             Normal Spectrum brginning with minimum :    0\n");
  23.   printf("                                                peak:    3\n");
  24.   printf("             When specifying -tri, the phases will be recalculated !\n");
  25.   printf("  -get file  get old spectra definition file (for reading\n");
  26.   printf("                              the spectrum phase and parity)\n");
  27.   printf("  -def file  specifies another definition file for output\n");
  28.   printf("  -v         print peak positions to stdout\n\n");
  29.   printf("the default output file is %s\n",TMPFILE);
  30.   printf("This file is read by the RVT routine and is organized as follows:\n");
  31.   printf("Each line consists of 4 integers.\n");
  32.   printf("The first specifies the Time 0 position.\n");
  33.   printf("The second can be +1 -1 +2 or -2 and specifies, if it is\n");
  34.   printf("     to be mirrored (-) and if it is a 0 degree (1)\n");
  35.   printf("     or 90 degree (2) spectrum\n");
  36.   printf("The last two numbers are used to specify a range for\n");
  37.   printf("the RVT routine and the exponential fit.\n");
  38.   printf("\n(C) Rainer Kowallik\n");
  39.   exit(0);
  40. }
  41.  
  42. main(argc,argv)
  43. int argc;
  44. char *argv[];
  45. {
  46. int xx,yy,n,m,i,spcmax,peakptr;
  47. int peaks[MAXPEAKS], phase[MAXPEAKS], xend[MAXPEAKS];
  48. char z[80],defnam[80];
  49. float x,y;
  50. FILE *fp;
  51.  
  52.    checkopt(argc,argv,"-dummy",z);
  53.  
  54.    spc= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  55.    err= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  56.    tim= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  57.  
  58.    strcpy(defnam,TMPFILE);
  59.  
  60.    guessphase(phase,tri);
  61.    getphase(defnam,phase,&tri,&range);
  62.  
  63.    if(checkopt(argc,argv,"-tri",z)) {
  64.       tri=atoi(z);
  65.       if(tri>3) {
  66.          printf("don't fool me , or I will fool you too !\n");
  67.          exit(0);
  68.       }
  69.       guessphase(phase,tri);
  70.    }
  71.    if(checkopt(argc,argv,"-def",z)) strcpy(defnam,z);
  72.  
  73.    if(checkopt(argc,argv,"-get",z)) getphase(z,phase,&tri,&range);
  74.  
  75.    spcmax=readspec(argv[1],spc,err,tim,z);
  76.  
  77.    peakptr = find_peaks(spc,peaks,MAXPEAKS,10,spcmax);
  78.    if(checkopt(argc,argv,"-v",z)) {
  79.       for(i = 0; i < peakptr; i++) {
  80.          n = peaks[i];
  81.          printf("peak at %5d height %f\n",n,spc[n]);
  82.       }
  83.    }
  84.    peakptr = calc_phases(spc,peaks,phase,peakptr,spcmax); /* !!!!!!!!!! */
  85.  
  86.    /* ---------------------------------------
  87.       finding end of usable data for spectra
  88.       --------------------------------------- */
  89.  
  90.    findusable(peaks,range,tri,phase,xend,peakptr);   
  91.    n=xend[0];
  92.    for(i=0;i<peakptr;i++) {
  93.       if(xend[i]<n) n=xend[i];
  94.    }
  95.    xend[peakptr-1] = n;                   /* minimum must be the last entry */
  96.  
  97.  
  98.    /* ---------------------------------------------------
  99.       writing data 
  100.    ------------------------------------------------------ */
  101.  
  102.    fp=fopen(defnam,"w");
  103.    for(i=0;i<peakptr;i++) {
  104.       fprintf(fp,"%d   %2d   %d   %d\n",peaks[i],phase[i],5,xend[i]);
  105.    }
  106.    fclose(fp);
  107.    printf("\ndata written to %s\n now start editor on this file !\n",TMPFILE);
  108.    free(spc); free(err); free(tim);
  109.    exit(0);
  110. }
  111.  
  112. iabs(x)
  113. int x;
  114. {
  115. if(x>=0) return(x);
  116. return(-x);
  117. }
  118.  
  119.  
  120. /* -----------------------------------------------
  121.    reading spectra definition file
  122.    ----------------------------------------------- */
  123. getphase(name,phase,tri,range)
  124. char name[];
  125. int phase[], *tri, *range;
  126. {
  127. int i,n,m,p[30];
  128. FILE *fp;
  129.  
  130.    fp=fopen(name,"r");
  131.    if(fp==NULL) {
  132.       return(0);
  133.    }
  134.    i=0; *tri=0;
  135.    while(!feof(fp)) {
  136.       fscanf(fp,"%d %d %d %d\n",&p[i],&phase[i],&n,&n);
  137.       i=i+1;
  138.    }
  139.    fclose(fp);
  140.    for(n=0;n<i;n++) if(phase[i]<0) *tri=2;   /* triangular spectrum ? */
  141.    if(phase[0]<0) *tri=1;        /* triangular, beginning with minimum */
  142.    n=p[0];
  143.    if(p[1]>(n+5)) {
  144.       n = iabs(n - p[1]);
  145.    } else {
  146.       n = iabs(n - p[2]);
  147.    }
  148.    i=1;
  149.    while(i<n) i= i << 1;
  150.    m = i >> 1;
  151.    if(iabs(n-m) < iabs(n-i)) {
  152.       n = m;
  153.    } else {
  154.       n = i;
  155.    }
  156.    *range = n;
  157. }
  158.  
  159. guessphase(phase,tri)
  160. int phase[], tri;
  161. {
  162. int i,n,m;
  163.  
  164.    n=1;
  165.    for(i=0;i<30;i++) {
  166.       n= -1 * n;
  167.       switch(tri) {
  168.       case 0:
  169.          if(n>0) phase[i]=1;
  170.          if(n<0) phase[i]=2;
  171.          break;
  172.       case 1:
  173.          if(n>0) {
  174.             phase[i++] = -1;
  175.             phase[i]=1;
  176.          }
  177.          if(n<0) {
  178.             phase[i++]= -2;
  179.             phase[i]=2;
  180.          }
  181.          break;
  182.       case 2:
  183.          if(n>0) {
  184.             phase[i++]=1;
  185.             phase[i]= -1;
  186.          }
  187.          if(n<0) {
  188.             phase[i++]=2;
  189.             phase[i]= -2;
  190.          }
  191.          break;
  192.       }
  193.    }
  194. }
  195.  
  196. /* ------------------------------------------------
  197.    finding usable data area
  198.    ------------------------------------------------ */
  199. findusable(peaks,range,tri,phase,xend,peakptr)
  200. int peaks[],phase[],xend[],
  201.     range,tri,peakptr;
  202. {
  203. int i,n,m,start,end;
  204.  
  205.    start = range / 2;
  206.    end = range;
  207.    if(tri>0) {
  208.       start = range / 4;
  209.       end = range / 2;
  210.    }
  211.    end = end - 3;
  212.    for(n=0;n<peakptr;n++) {  /* scanning spectra */
  213.       xend[n] = end;
  214.       if(phase[n] > 0) {    /* positive time axis */
  215.          m = peaks[n] + start;
  216.          for(i=start;i<end;i++) {
  217.             m = m + 1;
  218.             if(!outrange(spc[m],err[m],spc[m+1],err[m+1])) continue;
  219.             if(!outrange(spc[m],err[m],spc[m+2],err[m+2])) continue;
  220.             if(!outrange(spc[m-1],err[m-1],spc[m+1],err[m+1])) continue;
  221.             xend[n] = i - 2 ;
  222.          }
  223.       } else {
  224.          m = peaks[n] - start;
  225.          for(i=start;i<end;i++) {
  226.             m = m - 1;
  227.             if(!outrange(spc[m],err[m],spc[m-1],err[m-1])) continue;
  228.             if(!outrange(spc[m],err[m],spc[m-2],err[m-2])) continue;
  229.             if(!outrange(spc[m+1],err[m+1],spc[m-1],err[m-1])) continue;
  230.             xend[n] = i - 2 ;
  231.          }
  232.       }
  233.    }
  234. }
  235.  
  236. outrange(x1,e1,x2,e2)
  237. float x1,e1,x2,e2;
  238. {
  239. float a,b,c;
  240.  
  241.    a = (float) fabs((double)(x1-x2));
  242.    b = 3 * (e1 + e2);
  243.    if(a>b) return(TRUE);
  244.    return(FALSE);
  245. }
  246.  
  247. /* ------------------------------------------------------------------------
  248.             Here starts the Peak finder routine:
  249.    Algorithm:
  250.    1. calculate arithmetic mean of spectrum.
  251.    2. copy spectrum to work area.
  252.    3. Find Maximum of Work area.
  253.    4. find points left and right, which are below the arithmetic mean.
  254.    5. Set area between these points to 0
  255.    6. goto 3 until no points are above the mean value.
  256.    ------------------------------------------------------------------------ */
  257.  
  258. find_peaks(spc,peaks,maxpeaks,part,spcmax)
  259. float spc[];
  260. int peaks[];
  261. int maxpeaks;
  262. int part;
  263. int spcmax;
  264. {
  265. int peakptr, n, m, i;
  266. float *buffer, mean, low;
  267.  
  268.    mean = 0.0;
  269.  
  270.    buffer = (float *) calloc(spcmax+2,sizeof(float));
  271.    for(i = 0; i < spcmax; i++) {
  272.       mean += spc[i];
  273.       buffer[i] = spc[i];
  274.    }
  275.    mean = mean / spcmax;
  276.  
  277.    for(peakptr = 0; peakptr < maxpeaks; peakptr++) {
  278.       n = eliminate_max(buffer,mean,spcmax);
  279.       if(n < 0) break;
  280.       if(peakptr == 0) low = spc[n] / part;
  281.       if(spc[n] < low) {
  282.          peakptr--;
  283.          continue;
  284.       }
  285.       peaks[peakptr] = n;
  286.    }
  287.    sort(peaks,peakptr);
  288.    free(buffer);
  289.    return(peakptr);
  290. }
  291.  
  292. eliminate_max(spc,mean,spcmax)
  293. float spc[];
  294. float mean;
  295. int spcmax;
  296. {
  297. int i,n,m, x_max, lmin, rmin;
  298. float y_max;
  299.  
  300.    y_max = mean; x_max = -1;           /* find maximum */
  301.    for(i = 0; i < spcmax; i++) {
  302.       if(spc[i] > y_max) {
  303.          y_max = spc[i];
  304.          x_max = i;
  305.       }
  306.    }
  307.    if(x_max < 0) return(x_max);
  308.    lmin = x_max;                       /* find left minimum */
  309.    while(spc[lmin] > mean) lmin--;
  310.    rmin = x_max;                       /* find right minimum */
  311.    while(spc[rmin] > mean) rmin++;
  312.    for(i = lmin; i<rmin; i++) spc[i] = 0.0;  /* eliminate peak */
  313.  
  314.    return(x_max);
  315. }
  316.  
  317. sort(buffer,p)
  318. int buffer[];
  319. int p;
  320. {
  321. int i,n,m,z;
  322.  
  323.    for(n = 1; n < (p - 1); n++) {
  324.       i = 1;
  325.       for(m = 0; m < (p - n); m++) {
  326.          if(buffer[m] > buffer[m+1]) {
  327.             z = buffer[m];
  328.             buffer[m] = buffer[m+1];
  329.             buffer[m+1] = z; i = 0;
  330.          }
  331.       }
  332.       if(i == 1) break;
  333.    }
  334. }
  335.  
  336. calc_phases(spc,peaks,phase,peakptr,spcmax)
  337. float *spc;
  338. int *peaks;
  339. int *phase;
  340. int peakptr;
  341. int spcmax;
  342. {
  343. int i,n,m;
  344. float mean;
  345. int x1,x2;
  346. float y1,y2;
  347.  
  348.    mean = 0.0;
  349.    for(n = 0; n < spcmax; n++) mean += spc[n];
  350.    mean = mean / spcmax;
  351.  
  352.    for(n = 0; n < peakptr; n++) {
  353.       x1 = peaks[n] - 50;
  354.       x2 = peaks[n] + 50;
  355.       y1 = spc[x1];
  356.       y2 = spc[x2];
  357.       if((y1 > mean) && (y2 > mean)) {    /* triangular spectrum */
  358.          for(i = peakptr; i >= n; i--) peaks[i+1] = peaks[i];
  359.          peakptr = peakptr + 1;
  360.          phase[n] = -1;
  361.          n = n + 1;
  362.          phase[n] = 1;
  363.       }
  364.       if((y1 > mean) && (y2 < mean)) {    /* normal inverted spectrum */
  365.          phase[n] = -1;
  366.       }
  367.       if((y1 < mean) && (y2 > mean)) {    /* normal, NOT inverted spectrum */
  368.          phase[n] = 1;
  369.       }
  370.    }
  371.    return(peakptr);
  372. }
  373.  
  374.